
---
title: "Latta et al 2017 PeerJ Rarefactions"
author: "brouwern@gmail.com"
date: "May 3, 2017"
output: html_document
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


## Rarefaction analysis

General workflow:

* Load mist nest and point count data  
* Seperate Aug and Jan dataframes
* Turn dataframes into lists compatible with iNEXT
* Pool data across years
* Run iNEXT to get rarefaction curves for **each year** (not used in final MS)
* Run iNEXT to get rarefaction curves for pooled data
* Plot



## Load data

### Load processed mist net data


Set the working directory and load the data.

```{r, message=FALSE, warning=FALSE}
#setwd

my.dir <-  "C:/Users/lisanjie2/Dropbox/Latta_et_al_2017_PeerJ_FINAL"

setwd(my.dir)

#load data
mistnet <- read.csv("./DATA/Latta_et_al_mistnet_wide.csv")


dim(mistnet)               #200    41; wide data
```

### Load point count data

```{r}

#setwd
my.dir <-  "C:/Users/lisanjie2/Dropbox/Latta_et_al_2017_PeerJ_FINAL"

setwd(my.dir)

#load data
pc <- read.csv("./DATA/Latta_et_al_point_counts_long.csv")

```



## Load libraries
```{r, echo=F}
library(vegan)
library(iNEXT)
library(ggplot2)
library(reshape2)
library(doBy)
```


## Subset wide mist nest data by month

NB: Aug & Jan don't have exactly same years; one has 2 fewer years

```{r}
#names of columns
year.names <- c("n.05","n.06","n.07","n.08",         
"n.09","n.10","n.11","n.12",
"n.13","n.14")

#subset by year
CR.mistnet.yrXspp.Jan <-subset(mistnet, month == "Jan", select = year.names)
CR.mistnet.yrXspp.Aug <-subset(mistnet, month == "Aug", select = year.names)

```


### Create list of mist nest data for iNEXT - January

iNEXT works with LISTS, not dataframes.  This code takes a dataframe and makes a seperate list element for each year.

```{r}

#create dummy list to hold processed data
CR.mistnet.Jan.list <- as.list(1:dim(CR.mistnet.yrXspp.Jan)[2])

#Loop over each COLUMN and put data for each YEAR into seperate element of the list
for(i in 1:dim(CR.mistnet.yrXspp.Jan)[2]){

x.use <- CR.mistnet.yrXspp.Jan[,i]
x.use <- x.use[which(x.use !=0)]
CR.mistnet.Jan.list[[i]] <- x.use
}


#each list is named by the year
names(CR.mistnet.Jan.list) <- paste("year",c(2005:2014),sep = ".")

```

### Create list of mist nest data for iNEXT - August

August has fewer data points than January (2 fewer years)

```{r}
#remove years when data not collected
CR.mistnet.yrXspp.Aug2 <- CR.mistnet.yrXspp.Aug[,-c(2, 8)]

CR.mistnet.Aug.list <- as.list(1: dim(CR.mistnet.yrXspp.Aug2)[2])

for(i in 1:dim(CR.mistnet.yrXspp.Aug2)[2]){

x.use <- CR.mistnet.yrXspp.Aug2[,i]
x.use <- x.use[which(x.use !=0)]
CR.mistnet.Aug.list[[i]] <- x.use
}

names(CR.mistnet.Aug.list) <- paste("year",c(2005,2007:2011,2013,2014),sep = ".")

```

The lists look like this.  Again, each list element is a seperate year.  Species information ID not included or needed
```{r}
CR.mistnet.Aug.list[1:3]
```



# Pool all years of mist net data

This code takes the dataframe for each month and calcualtes the **total number of individuals** captured for each species across **all year** of the study.

This is what is used in the final MS.
```{r}
#make pooled column
## Aug
CR.mistnet.yrXspp.Aug$Aug.pool <- apply(CR.mistnet.yrXspp.Aug,1,sum, na.rm = T)

## Jan
CR.mistnet.yrXspp.Jan$Jan.pool <- apply(CR.mistnet.yrXspp.Jan,1,sum, na.rm = T)

```


Make a list with pooled data from each month
```{r}
# Make list
Jan.vs.Aug.pooled.list <- list(January = CR.mistnet.yrXspp.Jan$Jan.pool,
August =CR.mistnet.yrXspp.Aug$Aug.pool)

#the list looks like this
str(Jan.vs.Aug.pooled.list)
```






## iNEXT rarefacation of mist net data by year

###  iNEXT Janury mist net rarefaction by year
```{r}
CR.mistnet.iNEXT.Jan <- iNEXT(CR.mistnet.Jan.list, 
q=0, 
datatype="abundance"#, 
#endpoint=500
)

```

### iNEXT August mist net rarefaction by year
```{r}
CR.mistnet.iNEXT.Aug <- iNEXT(CR.mistnet.Aug.list, 
q=0, 
datatype="abundance"#, 
#endpoint=500
) 
```



##  iNEXT rarefaction on pooled mistnet data 

This does rarefaction of the pooled data.  This is what is used in the final MS.

```{r}
CR.mistnet.iNEXT.Jan.vs.Aug.pooled <- iNEXT(Jan.vs.Aug.pooled.list, 
q=0, 
datatype="abundance"#, 
#endpoint=500
)

```



## Subset point count data data by month and put into list

Again, iNEXT takes data only as a list, not dataframe
```{r}
pc.aug <- subset(x = pc,month == "aug")
pc.jan <- subset(x = pc,month == "jan")

pc.list<- list(August = pc.aug$count,
               January = pc.jan$count)
```



## Run iNEXT on point count data
```{r}
pc.iNEXT <- iNEXT(pc.list, 
             q=c(0), 
             datatype="abundance"#, 
             #endpoint=500
     )
```



## Plot point count Rarefaction Output
```{r}
ggiNEXT(pc.iNEXT, type=1,se = T) + 
  ylab("Species Richness") + 
  theme_bw() +
  ggtitle("iNEXT Rarefaction-Point Counts") +
  #xlim(c(0,1500)) +
  annotate("text",x = 275, y = 120, 
           label = "Error bands = +/- 1 SE")


```


## Combine mist nest and point count data


### Mist net data
```{r}
net.Jan.fin <- CR.mistnet.iNEXT.Jan.vs.Aug.pooled$iNextEst$January
net.Aug.fin <- CR.mistnet.iNEXT.Jan.vs.Aug.pooled$iNextEst$August
net.Jan.fin$Month <- "January"
net.Aug.fin$Month <- "August"

net.iNEXt.final <- rbind(net.Jan.fin,
                         net.Aug.fin)
net.iNEXt.final$Method <- "Mist net"

```



### Point count data

                        
```{r}
pc.Jan.fin <- pc.iNEXT$iNextEst$January
pc.Jan.fin$Month <- "January"

pc.Aug.fin <- pc.iNEXT$iNextEst$August
pc.Aug.fin$Month <- "August"
pc.iNEXt.final <- rbind(pc.Jan.fin,
                        pc.Aug.fin)

pc.iNEXt.final$Method <- "Point count"


```


## Combine mist net and point count data
```{r}
iNEXT.final <- rbind(pc.iNEXt.final,
                     net.iNEXt.final)

```

## Set up data for plotting Chao1 estimates

This produces a plot wiith the Chao1 estiamtes and CIs for the point 
count data and the pooled mist net data
In the iNEXT output, "AsyEst" contains the Chao1 data


### Maka seperate object for each for each set of hill numbers (chao1 estiamtes)
```{r}
pc.Asy <- pc.iNEXT$AsyEst
net.Asy <- CR.mistnet.iNEXT.Jan.vs.Aug.pooled$AsyEst

```



### Label the appropriaste methods

```{r}
pc.Asy$Method <- "Point Count"
net.Asy$Method <- "Mist Net"
```



### Make combined dataframe
```{r}
both.methods.Asy <- rbind(pc.Asy, net.Asy)

```

### Chagne "site" label to "month"
```{r}
names(both.methods.Asy)[1] <- "Month"

```


### Melt data into long format

```{r}
library(reshape2)
        
both.methods.Asy.long <- melt(data = both.methods.Asy,
                                      measure.vars = c("Observed","Estimator"))
        
i.obs <- which(both.methods.Asy.long$variable == "Observed")
        
#remove SE etc from observed data points
both.methods.Asy.long[i.obs, c("s.e.","LCL","UCL")] <- NA

```



        
        
## Plot the chao1 estiamtes/hill number fro the 2 methods next to each other

Not in MS

```{r}
pd <- position_dodge(0.25)
        ggplot(data = both.methods.Asy,
               aes(x = Method,
                   y = Estimator,
                   color = Month,
                   shape = Month)) +
          geom_point(  position = pd, 
                       size = 4) +
          geom_errorbar(aes(ymax = UCL,
                            ymin = LCL),
                        width = 0,
                        size = 1.15,
                        position = pd) +
          geom_point(data = both.methods.Asy,
                     aes(x= Method,
                         y = Observed),
                     position = pd,
                     size = 2,
                     color = "black") +
          
          facet_grid(facets = Diversity~.,scale = "free") +
          theme_bw() + 
          ylab("Species Diversity Estimate") +
          xlab("Sampling Method") +
          geom_vline(xintercept = 1.5) #+
        #annotate(geom = "text",x = 1.75)
```





## Rarefaction plots: Mist nest  vs Point counts

This figure appears in final MS

```{r}
i.interp <- which(iNEXT.final$method == "interpolated")
i.extrap <- which(iNEXT.final$method == "extrapolated")

iNEXT.final$Month2 <-iNEXT.final$Month

ggplot(data = iNEXT.final[i.interp,],
       aes(y = qD,
           x = m)) + 
  
  #Plot observed richness
  geom_point(data = iNEXT.final[-c(i.extrap,i.interp),
                                   ],
             aes(x = m,
                 y = qD,
                 color = Month2,
                 shape = Month2),
             size = 5) +
  
  #Plot "interpolated" data
  geom_line(aes(color = Month,
                linetype = Month),
            size = 2.25) +
  
  geom_ribbon(data = iNEXT.final[i.interp,],
              aes(x = m,
                  ymax = qD.UCL,
                  ymin = qD.LCL,
                  color = Month,
                  fill = Month),
              alpha = 0.3451,
              colour = 0) +

  #Plot extrapolated data
  geom_line(data =iNEXT.final[i.extrap,],
            aes( #linetype = Month2,
                group = Month2),
            size = 1.25)+
  geom_ribbon(data = iNEXT.final[i.extrap,],
              aes(x = m,
                  ymax = qD.UCL,
                  ymin = qD.LCL,
                  color = Month,
                  fill = Month),
              alpha = 0.3451,
              colour = 0)+
  
  
  
  #Misc
  facet_grid(facets = . ~ Method) +
  theme_bw() +
  xlab("") +
  ylab("") +
  xlim(0,2250) +
  
  #remove legend
  theme(legend.position = "none") + 
 
  #remove facet labeling and boxes around label
  theme(strip.background = element_blank(),
       strip.text.x = element_blank()) +
  
  #increase font size of axis labels
  theme(axis.text=element_text(size=18,face = "bold"))

```









## Make  Table of POOLED Data

### Point count data

```{r}

pc.iNEXt.AsyEst.final <- pc.iNEXT$AsyEst

pc.iNEXt.AsyEst.final$Method <- "Point count"

```


### Pooled Mist net data
```{r}
net.iNEXt.AsyEst.final <- CR.mistnet.iNEXT.Jan.vs.Aug.pooled$AsyEst

net.iNEXt.AsyEst.final$Method <- "Mist net"


iNEXT.AsyEst.final <- rbind(pc.iNEXt.AsyEst.final,
                            net.iNEXt.AsyEst.final)

```

Save
```{r}
#write.csv(iNEXT.AsyEst.final, 
          file = "COSTA_RICA_mist_net_and_PC_diversity_metrics_TABLE_vs1_Nov11.csv")
```

